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Abstract — For the non-Newtonian mantle rheology case the 2D thermal viscous dissipation-driven thermal convection in 
the mantle wedge above the Apulian lithospheric microplate subducting under the Euro-Asian plate is modeled numerically. 
The effects of the 410 km and 660 km phase transitions are taken into account. Within the framework of the model 
constructed the horizontal extent of the 2D heat flux anomaly observed in the rear of the Dinarides mountain belt 
corresponds to subduction velocity of -10 mm per year which is close to that observed with the help of geodetic means. In 
the case of non-Newtonian rheology the upwelling convective flow transporting heat to the Earth’s surface locates at the 
distance from the trench corresponding to the actually observed 2D heat flux anomaly, the velocity in the convective vortices 
being from -10 mm per year to -10 m per year for the water content in the mantle wedge from 0.3x10 to 3x10 weight 
percent respectively. The convection cell dimension is of the order of the horizontal scale of the heat flux anomaly observed 
in the Pannonia petroleum basin and Vardar zone. Upwelling mantle wedge convective flow is indicated to be able to 
provide the mantle wedge hydrocarbons transport to the Earth’s surface for the mantle wedge mantle content over 0.3xl0T 1 
weight percent. 
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I. Introduction 

According to [3], the subduction of the Apulian lithospheric microplate under the Dinarides, Pannonia basin and Vardar zone 
is sufficiently flat and during the last ~ 45 Ma occurred at the angle of ~ 25°, which remained unchanged during this time 
interval. The genesis of the Dinarides mountain belt (with the transversal horizontal extent of ~ 300 km) apparently is of the 
thrust and fold nature, associated with the former collision and subduction of the oceanic branches of the Neo-Tethis and 
Alpine-Tethis as the result of thrusting of the African plate under the Eastern and Western Europe during the last 55 - 35 Ma 

[3] . In [1] numerous papers are referred, containing contradictory estimates of the relative motions of the Apulian 
lithospheric microplate and Euro-Asian plates, made on the basis of seismic, geophysical and geodetic data. In fig. 3 in [Op. 
cit.] the velocity of subduction of the Apulian lithospheric microplate under the Euro-Asian one is seen to amount to ~5 - 8 
mmxa" 1 according to referred works, while in [1] this velocity is estimated to be of the order of ~ 5 mmxa according to 
calculations based on geodetic observations. In [12] the Pannonia petroleum basin and Vardar zone are noted to be the zones 
of the Middle Miocene extension occurred -14-11.6 Ma ago, which led to the lithosphere thinning, these zones being the 
back-arc basin characterized by the back-arc spreading. At that time the single mountain belt parallel to Apulian shore was 
split into Carpathians and Dinarides and the shallow Pannonian sea was formed, which existed approximately to 600 
thousand years ago. Presently the Pannonian oil- and gas-bearing basin is situated in this region. Here the conditions are 
clarified under which the centre of the back-arc spreading initiates as the result of convective instability, driven by the 
dissipative heat release in the mantle wedge above the subducting Apulian lithospheric microplate. 

According to [4, 8, 9] two types of dissipation-driven small-scale thermal convection in the mantle wedge are possible, viz. 
the 3D finger-like convective jets, raising to volcanic chain, and 2D transversal Karig vortices, aligned perpendicularly to 
subduction. These two types of convection are shown to be spatially separated due to the pressure and temperature 
dependence of mantle effective viscosity, the Karig vortices, if any of them formed, being located behind the volcanic arc 

[4] . There are contradictory judgments on the velocity of subduction of the Apulian lithospheric microplate under the Euro- 
Asian one, although the order of magnitude of the present-day subduction velocity (~ 10 mmxa -1 ) can apparently be regarded 
as established sufficiently reliably. The mountainous massif Dinarides locates parallel to the north-eastern shore of the 
Adriatic sea, and probably is of the thrust-and-fold nature. The 2D maximum of the heat flux anomaly of -100 mWxm 2 




International Journal of Engineering Research & Science (IJOER) 


ISSN, [2395-6992] 


[Vol-5, Issue-10, October- 2019] 


observed in the rear of the Dinarides massif in the Pannonia basin and the Vardar zone [10] can be assumed to owe its origin 
to the convective heat supply from the mantle wedge. Numerical modeling of the 2D convection, occurring in the mantle 
wedge in the form of the Karig vortices and presumably transporting heat upwards, allows to judge about the mean water 
content in the mantle wedge and to assume the mantle hydrocarbons to be transported to the Earth’s surface by the upwelling 
convective flows. Numerical convection models accounting for the pressure-, temperature- and stress-dependence of 
viscosity fit best to observational data in the case of non-Newtonian rheology at the mantle water content of ~ (0.3 - 3)1 10 1 
weight % for the velocity of subduction of ~ 10 mmxa 1 during the Middle Miocene. In [15] such rather a high water content 
(and even greater one, up to 3 weight %) can be observed in the mantle wedge in the mantle transition zone. The Middle 
Miocene subduction velocity of ~10 mmxa during the formation of the Pannonia basin is of the order of the observed 
presently, or, can somewhat exceed it because of the gradual diminution of the velocity of convergence of African and Euro- 
Asian plates. 

II. The Computer Thermo-Mechanical Model 

Thermo-mechanical model of the mantle wedge between the base of the overlying Euro-Asian plate and the upper surface of 
the Apulian lithospheric microplate subducting under the Euro-Asian one with a velocity V at an angle p is obtained for the 
infinite Prandtl number fluid as a solution of non-dimensional 2D hydrodynamic equations in the Boussinesq approximation 

(S 2 zz -d 2 xx )x ij x (5 2 zz -5 2 xx )x y/ +4 *d 2 xz rj x^y / =R a xT x -Ra^T^ m -Ra^xT^ (1) 

d l T=AT-(i// z x T x ) + ( y/ x x T z ) + (Di / Ra ) x (r, A , / 2x p + Q (2) 


For the stream-function yr and temperature T. Here ?j is dynamic viscosity, d and indices denote partial derivatives with 
respect to coordinates x (horizontal), z (vertical) and time t, A is the Laplace operator, T t4ln) and r (660) are volume ratios of the 
heavy phase at the 410 km and 660 km phase boundaries, the velocity components V x and V z are expressed through (/x as 


V x = y/ z ,V z = - y/ x (3) 

while non-dimensional Rayleigh number Ra, phase numbers Ra t4l> ", Ra (660) and dissipative number Di are 

Ra = [(a xp xgxd 3 x7j) / ( tj c x /)] = 5.55xl0 8 ; 

Ra (410 4(5p (410, x gxd 3 )/(ii c x /)]= 6.60xl0 8 ; (4) 

Ra (6m = [(5p (660) x gxd 3 ) / (Tjc x /)] = 8.50x10 s ; 

Di = |(a x gx d) / c p \ = 0.165, 

where a= 3xl0’ 5 K' 1 is the thermal expansion coefficient, p = 3.3 g x cm' 3 is the density, g is gravity acceleration, c p = 
1.2xl0 3 J x kg'xK’ 1 is specific heat capacity at constant pressure, 7) = 1950°K is the temperature at the base of mantle 
transition zone (MTZ) at depth 660 km regarded the lower boundary of the model domain, Q = 6.25xl0' 4 mWx m 3 is the 
volumetric heat generation in the crust, r ik is the viscous stress tensor, d = 660 km is vertical dimension of the modeled 
domain, rj s = 10 18 Paxs is the viscosity scaling factor, /= 1 mm 2 xs~* is thermal diffusivity, Sp (4U)> = 0.07pand Sp (66i>) = 
0.09 p are the density changes at the 410 km and 660 km phase boundaries respectively. In (1), (2) the scaling factors for time 
t, coordinates x and z, stresses r ik and the stream-function y/ are (d 2 x % _1 ), d, (tj s x % x d' 2 ) and % respectively. Assuming 
rheology be linear for the diffusion creep deformation mechanism dominating in the mantle at depths over - 200 km [2], we 
accept the temperature- and lithostatic pressure p dependent viscosity as [15]: 

il = (p 12 x A) x ( h/b*) m x { exp [ (E* + pxV*) i(RxT)]} (5) 

where for “wet” olivine A = 5.3xl0 15 s" 1 , m = 2.5, the grain size h = 10" 1 - 10 mm, b = 5xl0' 8 cm is the Burgers vector [14], 
E f = 240 kJxmol" 1 is activation energy, U =5xl0 3 mm 3 'mol' 1 is activation volume, p = 300 GPa is the shear modulus 
normalizing factor, R is the gas constant. At the constants chosen and the grain size h =1.6 mm, non-dimensional viscosity 
also denoted t] is 

i] = 5x10" 7 x exp {[14.8 + 6.72 x(l-z)]/7’} (6) 

where T is non-dimensional temperature, non-dimensional z normalized by d is pointing upwards from the MTZ base and x is 
pointing against subduction along the MTZ base. The aspect ratio of the model domain is 1:2.25 thus the subduction angle 
being P = 24° if subduction is assumed to take place along the model domain diagonal. Non-dimensional subduction velocity 
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V = 10 mmxfl 1 normalized by (% x d~ ') equals V = 0.208xl0 3 , i.e. non-dimensional velocity components of subducting 
Adriatic micro-plate are V x = - 0.190xl0 3 and Vz = — 0.085xl0 3 . 

To check as to how the estimate of velocity of subduction of the Adriatic micro-plate is sensitive to the accepted linear 
rheological law here we make extra computations for the non-Newtonian rheology, in which case the viscosity formulae (5)- 
(6) are rewritten as: 

ij = (1 / 2 x A x C w r x r"-‘) x (h / b*) m x { exp [ (E* + pxV*) / (R xT) ] } (7) 

where according to [13] for “wet” olivine n = 3, r = 1.2, m = 0, x = (t 2 ,*-) 1/2 , E*= 480 kJ'mol" 1 , V* = llxlO 3 mm 3 xmol' 1 , A = 
10 2 c 'xJMPa)’", C„, > 10" 3 for “wet” olivine is the weight water concentration (in %%). It should be noted the constants in (7) 
vary considerably in the papers referred to by [13] and heretofore we gave averaged values of constants. At C w =10" 3 on 
accounting for 

Tik = (4 X ?7 2 ) x [ (y/ zz - y/ xx f I 2 + 2 x y/J ] (8) 

non-dimensional viscosity is 

? 1 = { 1.0 / [ ( ¥j7 - yjt 2 + 2 x y/J ] 1/3 } xexp{ [10.0 + 5.0 x (1 - z)] / T } (9) 

Following [13] we assume the phase functions r (,) as 

r (,) = (l/2)x{l -th [z-z (l> (T)]/w (l> };z il) (T) = z 0 (l> - {[fx (T - T o 0> )] / (p x g)} (10) 

where the signs are changed as z-axis is pointing upwards, z <l> (T) is the depth of the /-th phase transition (/ = 410, 660), z 0 (l> 
and Tj 11 are the averaged depth and temperature of the /-th phase transition, y <410) = 3 MPaxK and y <660) = -3 MPaxK 1 are 
the slopes of the phase equilibrium curves, w {U is the characteristic thickness of the /-th phase transition, 7]/ 410 ’ =1800° K, 
7’„ <66 o , = 1950 0 K are the mean phase transition temperatures. The heats of phase transitions are neglected in (2) as insignificant 
in the case of developed convection as in [13]. From (10) it follows 

T x (,) = - (y (,) / 2 xpxgxw (l> ) x T x x ch~ 2 { [(z - zj n + f x (T-T 0 ,l> )) / (pxg)] / w ,n ] (11) 

Where from it is clear the phase transition with y (l> > 0 facilitates convection (at / = 410), while the phase transition with y ,,J < 
0 hinders convection (at / = 660). In non-dimensional form zj 4m = 0.38, z o (660) = 0, w (l> = 0.05, y <410) = 2.5xl0 9 , y (660) = - 
2.5xl0 9 , T„ (410) = 0.92, t o (660) = 1 and in (1) 

r r (,) =-(8'p w y f, V2xp Ra w xw (I )xT x xch- 2 {[z-zJ l> W\Sp a) /pRad)x(T-T 0 (l> )] / w {,) } (12) 

Equations (l)-(2) are solved for the isothermal horizontal and vertical boundaries regarded no-slip impenetrable ones except 
for the “windows” for in- and outgoing subducting plate, where the plate velocity is specified. Vertical boundary distant from 
subduction zone is assumed penetrable at right angle; the latter boundary condition appears not too imposing in the case of 
very flat subduction. Q in (2) is non-zero in the continental and oceanic crust 40 and 7 km thick. Initial vertical boundaries 
temperature is calculated for the half-space cooling model for 10 9 yr and 10 8 yr for Euro-Asian (continental) plate and 
Apulian lithospheric microplate (paleooceanic plate) respectively. It is convenient to express dimensionless T 1K in (2) 
through the stream-function (/xas in (8): 

Tik 2 = (4xf)x[(y zz - y/ xx ) 2 / 2 + 2 x y / xz 2 ] 

III. Results and Discussion 

Assuming the heat flux maximum q is formed above the convective flow, upwelling towards the Pannonia basin and Vardar 
zone, and the convection cell dimension is equal to horizontal scale of the heat flux anomaly zone, the convection cell 
dimension can be estimated of ~ 300 km. To more accurately compute the consistent model of small-scale convection in the 
mantle wedge between the overriding Euro-Asian plate and subducting Apulian lithospheric microplate it is necessary from 
the computational point of view first to specify in (l)-(2) vanishing non-dimensional numbers /to—>0, Di = 0, i.e. to ignore 
convection and viscous dissipation. This approach is applied as convection with Rci and Di (4) passes through rather vigorous 
stages, and the time steps in integrating (1)—(2) become too small thus making it difficult to model the thermal structure of 
the plates. Solving (l)-(2) by the finite element method in space on the grid 104x104 and the 3-rd order Runge-Kutta method 
in time one obtains for Ra— >0, Di = 0 and l'= 10 mm a year non-dimensional quasi steady-state i// and T shown in fig. 1, 
where the streamlines are depicted with the step 5 and the isotherms with an interval 0.05. 
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FIGURE 1. Schematic cross section of the region of subduction of the Apulian lithospheric microplate under 
the Euro-Asian plate with no effects of viscous dissipation and convection taken into account. (1) - the 
quasi steady-state distribution of non-dimensional temperature and (2) - the quasi steady-state distribution 
of stream-function, the streamlines above the subducting slab correspond to the mantle wedge flows (“A” 

and “B”) induced by subduction. 


Subducting plate was considered rigid, while the viscosity at the zone of plates friction (at temperatures below 1200° K) was 
reduced by 2 orders of magnitude as compared to (7). The latter viscosity reduction at the plates contact zone accounts for 
lubrication effected by deposits partially entrained by the subducting plate. Such a lubrication prevents the overriding Euro- 
Asian plate from gluing to the subducting one [9]. Fig. 1 shows the results of computation for the formulae (7) - (9) for non- 
Newtonian rheology case for the water content C w = 3x10"* weight %. The velocity V= 10 mm per year is chosen as resulting 
in the best fit of the model convective zone size to horizontal extent of the observed heat flux anomaly in Pannonia basin and 
the Vardar zone. The Apulian lithospheric microplate subducting with a given velocity V is considered rigid and is shown in 
Fig. 1(2) by the equidistant diagonal streamlines. The induced mantle wedge flow above the subducting plate is seen to occur 
in the form of two vortices A and B (located one above another), the latter 2 vortices being considerably compressed in the 
vertical direction and the upper one (with y/ > 0) revolves clockwise while the lower one (with i// < 0) revolves 
counterclockwise. In Fig. 1(2) the upper induced flow “A” is seen to be firmly pressed to the subducting Apulian lithospheric 
microplate, the strain rate in the zone of contact of the opposite flows (i.e. flow “A” and subducting slab) being very high 
thus resulting in the viscosity (7) drop by several orders of magnitude. It should be noted that in the case of non-Newtonian 
rheology the greater the subduction angle the more extensive is this contact zone, in which the dissipative heat release mainly 
occurs. This may serve the reason why the Karig vortices (and the resulting back-arc spreading) are formed in the zones of 
comparatively steep subduction. The opposite flow “A” in Fig. 1(2) apparently is induced by the flow “B”, forced by 
subducting plate. 

Assuming Ra = 5.55xl0 8 and Di = 0.165, i.e. switching dissipation and convection on, taking into account the effects of 
phase transitions, from (1)—(2) the convection at C w = 3x10 1 weight % is found to destroy the induced mantle flows in the 
mantle wedge during the time interval ~ 0.6xl0' 6 (in dimensional form ~ 0.1 Ma) and to assume the quasi steady-state form 
shown in Fig. 2, in which the streamlines in the convection vortices are depicted with the interval 4xl0 4 . 
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FIGURE 2. Quasi steady-state non-dimensional stream-function distribution in the zone of subduction of the 
Apulian lithospheric microplate under the Euro-Asian plate with the effects of dissipative heating and 
convection taken into account for non-Newtonian rheology for the water content C w =3xl0 _1 weight %. 
Parallel equidistant streamlines represent subducting Apulian lithospheric microplate. Arrow “A” shows 
possible direction of the upward convective transport of dissipative heat and mantle hydrocarbons to the 

Earth’s surface. 

These convective vortices are seen actually to correspond to a single convection cell aroused at subduction velocity V = 10 
mmxa '. The latter convection cell dimension is of the order of ~ 300 km, i.e. is very close to the observed horizontal extent 
of the heat flux anomaly observed in the Pannonia basin and Vardar zone. The velocity in convective vortices may exceed 10 
mxfif 1 for the water content of 3 x 10 weight %. The direction of a possible upward transport of mantle hydrocarbons and 
dissipative heat to the Earth’s surface is shown by the arrow “A”. It should be noted that in the case of Newtonian rheology 
the convection in the mantle wedge cannot be aroused at the velocity of subduction of 10 mmxT 1 and subduction angle of 
25°. 

In Fig. 3 the quasi steady-state stream-function y /is shown for the mantle water content of 0.3x10"' weight %, in which case 
the dissipation-driven convection is aroused in essentially a single vortex with a characteristic velocity of -10 mmxa" 1 . 

Distance, km 



FIGURE 3. Quasi steady-state non-dimensional stream-function distribution in the zone of subduction of the 
Apulian lithospheric microplate under the Euro-Asian plate with the effects of dissipative heating and 
convection taken into account for non-Newtonian rheology for the water content C w = 0.3x10 ' weight %. 
Parallel equidistant streamlines represent subducting Apulian lithospheric microplate. Arrow “A” shows 
possible direction of the upward convective transport of dissipative heat and mantle hydrocarbons to the 

Earth’s surface. 

This convective flow may serve a means of the upward mantle hydrocarbons and dissipative heat transport, which direction 
is shown in Fig. 3 by the arrow “A”. Such a convective flow may cause the back-arc spreading in the Pannonia basin and the 
Vardar zone. At the mantle water content C„, < 0.3x10 1 weight percent the dissipation-driven convection cannot be aroused 
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in the mantle wedge even in the case of the non-Newtonian rheology for subduction velocity of ~10 mmxa 1 and angle of 
25°. 

It is worth noting that in the case of Newtonian rheology the 2D transversal convective rolls in the mantle wedge, as in Fig. 
2, can be formed only at sufficiently small angles of subduction. Thus, at /? = 30° the transversal convective rolls are not 
formed even at the velocity of subduction of 100 mmxrT 1 [5, 9]. In the case of the non-Newtonian rheology the transversal 
rolls (2D Karig vortices) can be aroused at greater subduction angles and sufficiently small subduction velocities owing to 
viscous friction in the zone of contact of the opposite induced flow (“A” in Fig. 1) and subducting slab. It should be noted 
that numerous thermo-mechanical mantle models in the zones of subduction (see, e.g. [8, 9] and the vast number of 
references there) showed convection in the form of transversal rolls never to occurred as the models with extremely small 
subduction angle and sufficiently great subduction velocity were not investigated. 

IV. Conclusion 

The size of the cell of 2D mantle wedge dissipation-driven convection in the case of the realistic non-Newtonian rheology 
equals ~ 300 km at the subduction velocity 10 mmxa -1 , in which case a single convection cell is aroused. This explains the 
formation and horizontal extent of the only 2D heat flux anomaly observed in the rear of the Dinarides. The water content 
sufficient for the 2D convection to be aroused is ~ 3x10"* weight %, or, alternatively, it is -3x]0 2 weight %, but the 2D 
convection is aroused as a single Karig vortex. The velocity in convective vortices in the non-Newtonian rheology case is 
~10 m per year at the water content C w ~ 3x10 weight percent and □ 10 mm per year at the water content C w ~ 0.3x10 
weight percent in the mantle wedge. The upwelling convective flow may be sufficient to provide upward transport of mantle 
wedge hydrocarbons to the Earth’s surface. 
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